Fit tweedie models to biomass density and extract AIC

Author

Max Lindmark

Published

2024-02-19

Intro

Fit Tweedie models to biomass density of cod, flounder, plaice and dab (juveniles and adults) between 1993-2020 in the Baltic Sea using sdmMTB fit different oxygen and temperature-related covariates, compare AIC. Select the “best” covariate for further trend and velocity analysis.

Load packages & source functions

# Load libraries, install if needed
pkgs <- c("tidyverse", "readxl", "tidylog", "RCurl", "devtools", "patchwork",
          "kableExtra", "viridis", "RColorBrewer", "here", "sdmTMBextra") 

if(length(setdiff(pkgs,rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library, character.only = T))

# Packages not on CRAN or dev version
# remotes::install_github("pbs-assess/sdmTMB", dependencies = TRUE)
library(sdmTMB)

# Source code for map plots
# You need: # devtools::install_github("seananderson/ggsidekick") # not on CRAN; library(ggsidekick)
devtools::source_url("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/R/functions/map-plot.R")
options(ggplot2.continuous.colour = "viridis")

# Set path
home <- here::here()

Read data

# Read data
d <- #readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/catch_clean.csv") %>% %>% 
  readr::read_csv(paste0(home, "/data/clean/catch_clean.csv")) %>%
  rename(X = x, Y = y) %>% 
  pivot_longer(c(cod_adult, cod_juvenile, dab_adult, dab_juvenile, flounder_adult,
                 flounder_juvenile, plaice_adult, plaice_juvenile),
               names_to = "group", values_to = "density") %>% 
  separate(group, into = c("species", "life_stage"), remove = FALSE) %>% 
  drop_na(depth, temp, oxy, sal, density)
Rows: 12385 Columns: 25
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): country, haul_id, ices_rect, month_year
dbl (21): year, lat, lon, quarter, month, sub_div, oxy, temp, sal, depth, x,...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
rename: renamed 2 variables (X, Y)

pivot_longer: reorganized (cod_adult, cod_juvenile, dab_adult, dab_juvenile, plaice_adult, …) into (group, density) [was 12385x25, now 99080x19]

drop_na: removed 9,570 rows (10%), 89,510 rows remaining
# Read metabolic parameter estimates and left_join
mi_pars <- #readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/mi_params.csv") %>%
  readr::read_csv(paste0(home, "/data/clean/mi_params.csv")) %>% 
  dplyr::select(n_o2, E_o2, A0_o2, species) #  TODO: remove the extra columns for pressure based parameters
Rows: 4 Columns: 10
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (1): species
dbl (9): temp, po2, o2, A0_o2, A0_po2, n_po2, E_po2, n_o2, E_o2

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# TODO: for now we'll use plaice parameters for flounder, see "00_estimate_mi_params.Rmd"
mi_pars
# A tibble: 4 × 4
    n_o2  E_o2  A0_o2 species 
   <dbl> <dbl>  <dbl> <chr>   
1 -0.300 0.294 0.0451 plaice  
2 -0.300 0.294 0.0249 cod     
3 -0.300 0.294 0.0371 dab     
4 -0.300 0.294 0.0120 flounder
mi_pars <- mi_pars %>% 
  mutate(A0_o2 = ifelse(species == "flounder",
                        filter(mi_pars, species == "plaice")$A0_o2,
                        A0_o2))
filter: removed 3 rows (75%), one row remaining
mutate: changed one value (25%) of 'A0_o2' (0 new NA)
d <- left_join(d, mi_pars, by = "species")
left_join: added 3 columns (n_o2, E_o2, A0_o2)
           > rows only in x        0
           > rows only in y  (     0)
           > matched rows     89,510
           >                 ========
           > rows total       89,510
# Read size csv to calculate the metabolic index
sizes <- #readr::read_csv("https://raw.githubusercontent.com/maxlindmark/spatial-metabolic-index/main/data/clean/sizes.csv") %>% 
  readr::read_csv(paste0(home, "/data/clean/sizes.csv")) %>%
  mutate(group = paste(species, name, sep = "_")) %>% 
  dplyr::select(group, B)
Rows: 8 Columns: 5
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (2): species, name
dbl (3): mat_w, max_w, B

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
mutate: new variable 'group' (character) with 8 unique values and 0% NA
d <- left_join(d, sizes, by = "group")
left_join: added one column (B)
           > rows only in x        0
           > rows only in y  (     0)
           > matched rows     89,510
           >                 ========
           > rows total       89,510
# Drop dab!
# d %>%
#   filter(group == "dab_juvenile" & quarter == 1) %>%
#   mutate(pres = ifelse(density > 0, "1", "0")) %>%
#   ggplot(aes(X, Y, color = pres)) +
#   geom_point(size = 0.3) +
#   coord_fixed() +
#   facet_wrap(~year)
# 
# d %>%
#   filter(group == "dab_juvenile" & quarter == 4) %>%
#   mutate(pres = ifelse(density > 0, "1", "0")) %>%
#   ggplot(aes(X, Y, color = pres)) +
#   geom_point(size = 0.3) +
#   coord_fixed() +
#   facet_wrap(~year)

d <- d %>% filter(!species == "dab")
filter: removed 23,124 rows (26%), 66,386 rows remaining

Calculate the metabolic index

# Oxygen is ml/L, We want micro mol/L. 1 ml/l = 10^3/22.391 = 44.661 micro mol/l
d$oxy_si <- (d$oxy * (10^3)) / 22.391

# Calculate metabolic indices for a given mass and the temperature and oxygen in data
# Line 123 in https://github.com/fate-spatialindicators/SDM_O2/blob/master/code/mi_functions.R
kb <- 0.000086173324 # Boltzmann's constant
t_ref <- 15 # arbitrary reference temperature

# Calculate the metabolic index
d <- d %>%
  mutate(inv_temp = (1/(temp + 273.15) - 1/(t_ref + 273.15)),
         phi = A0_o2*(B^n_o2)*oxy_si * exp((E_o2/kb)*inv_temp)) %>% 
  drop_na(phi)
mutate: new variable 'inv_temp' (double) with 12,194 unique values and 0% NA
        new variable 'phi' (double) with 65,602 unique values and 0% NA
drop_na: no rows removed

Scale variables

d <- d %>% 
  group_by(group) %>% 
  mutate(phi_sc = as.vector(scale(phi)),
         temp_sc = as.vector(scale(temp)),
         temp_sq = temp_sc^2,
         oxy_sc = as.vector(scale(oxy)),
         oxy_sq = oxy_sc^2,
         sal_sc = as.vector(scale(sal)),
         depth_sc = as.vector(scale(depth)),
         depth_sq = depth_sc*depth_sc) %>%
  mutate(quarter_f = as.factor(quarter),
         year_f = as.factor(year)) %>% 
  ungroup()
group_by: one grouping variable (group)
mutate (grouped): new variable 'phi_sc' (double) with 65,602 unique values and 0% NA
                  new variable 'temp_sc' (double) with 65,576 unique values and 0% NA
                  new variable 'temp_sq' (double) with 65,576 unique values and 0% NA
                  new variable 'oxy_sc' (double) with 65,568 unique values and 0% NA
                  new variable 'oxy_sq' (double) with 65,568 unique values and 0% NA
                  new variable 'sal_sc' (double) with 65,586 unique values and 0% NA
                  new variable 'depth_sc' (double) with 25,364 unique values and 0% NA
                  new variable 'depth_sq' (double) with 25,364 unique values and 0% NA
mutate (grouped): new variable 'quarter_f' (factor) with 2 unique values and 0% NA
                  new variable 'year_f' (factor) with 29 unique values and 0% NA
ungroup: no grouping variables
# Seems like a non-linear relationship with temperature, but a quadratic largely does the job??
ggplot(d, aes(temp_sc, density)) +
  geom_smooth(method = "lm", color = "grey50", se = FALSE) +
  geom_smooth(method = "gam", formula = y~s(x, k=4), color = "steelblue3", se = FALSE) +
  geom_smooth(method = "lm", formula = y ~ x + I(x^2), color = "tomato3", se = FALSE) +
  facet_wrap(~group, scales = "free")
`geom_smooth()` using formula = 'y ~ x'

ggplot(d, aes(oxy_sc, density)) +
  geom_smooth(method = "lm", color = "grey50", se = FALSE) +
  geom_smooth(method = "gam", formula = y~s(x, k=4), color = "steelblue3", se = FALSE) +
  geom_smooth(method = "lm", formula = y ~ x + I(x^2), color = "tomato3", se = FALSE) +
  facet_wrap(~group, scales = "free")
`geom_smooth()` using formula = 'y ~ x'

ggplot(d, aes(phi_sc, log(density + 1))) +
  geom_smooth(method = "lm", color = "grey50", se = FALSE) +
  geom_smooth(method = "gam", formula = y~s(x, k=4), color = "steelblue3", se = FALSE) +
  geom_smooth(method = "lm", formula = y ~ x + I(x^2), color = "tomato3", se = FALSE) +
  facet_wrap(~group, scales = "free")
`geom_smooth()` using formula = 'y ~ x'

# Plot mi by species and life stage
pal <- brewer.pal(name = "Dark2", n = 8)[c(7, 5, 3)]

d %>% 
  mutate(species = str_to_title(species),
         life_stage = str_to_title(life_stage)) %>% 
  ggplot(aes(phi, fill = species, color = species)) + 
  geom_density(alpha = 0.2) + 
  geom_vline(xintercept = 0, linetype = 2) + 
  coord_cartesian(expand = 0) + 
  facet_wrap(~ life_stage, ncol = 1) +
  theme_sleek(base_size = 9) + 
  theme(legend.position = c(0.85, 0.85)) +
  labs(x = "Metabolic index (\u03C6)", y = "Density", fill = "Species", color = "Species") + 
  scale_color_manual(values = rev(pal), name = "") +
  scale_fill_manual(values = rev(pal), name = "")
mutate: changed 66,386 values (100%) of 'species' (0 new NA)
        changed 66,386 values (100%) of 'life_stage' (0 new NA)

ggsave(paste0(home, "/figures/supp/mi_histogram.pdf"), width = 11, height = 11, units = "cm", device = cairo_pdf)

Fit models and save AIC

Fit models!

data_list_aic <- list()

for(i in unique(d$group)) {  
    
    dd <- d %>% filter(group == i)
    
    mesh <- make_mesh(dd, xy_cols = c("X", "Y"), cutoff = 15) # Need specific meshes because not for all hauls did we record all species
    
    print(
    ggplot() +
      inlabru::gg(mesh$mesh) +
      coord_fixed() +
      geom_point(aes(X, Y), data = dd, alpha = 0.2, size = 0.5) +
      annotate("text", -Inf, Inf, label = paste("n knots = ", mesh$mesh$n), hjust = -0.1, vjust = 2) + 
      labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(i, "_", " ")))
    )
    
    print(i)
    
    # 0. Linear oxygen and temperature, spatially varying quarter
    m0 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + temp_sc + oxy_sc,
                 data = dd,
                 mesh = mesh,
                 family = tweedie(link = "log"),
                 spatiotemporal = "IID",
                 # this setting for spatial field(s) works better for majority of groups compare to single sf
                 spatial = "off",
                 spatial_varying = ~0 + quarter_f, # if enabled, keep quarter_f in main formula; 
                 time = "year")
    print("m0")
    sanity(m0)
    dd$m0_res <- residuals(m0)
    
    # 1. Linear oxygen and squared temperature, spatially varying quarter
    m1 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + temp_sc + temp_sq + oxy_sc,
                 data = dd,
                 mesh = mesh,
                 family = tweedie(link = "log"),
                 spatiotemporal = "IID",
                 spatial = "off",
                 spatial_varying = ~0 + quarter_f,
                 time = "year")
    print("m1")
    sanity(m1)
    dd$m1_res <- residuals(m1)
    
    # 2. interaction between linear oxygen and temperature, drop quadratic
    m2 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + temp_sc * oxy_sc,
                  data = dd,
                  mesh = mesh,
                  family = tweedie(link = "log"),
                  spatiotemporal = "IID",
                  spatial = "off",
                  spatial_varying = ~0 + quarter_f,
                  time = "year")
    print("m2")
    sanity(m2)
    dd$m2_res <- residuals(m2)

    # 3. breakpoint oxygen
    m3 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + temp_sc + temp_sq + breakpt(oxy_sc),
                 data = dd,
                 mesh = mesh,
                 family = tweedie(link = "log"),
                 spatiotemporal = "IID",
                 spatial = "off",
                 spatial_varying = ~0 + quarter_f,
                 time = "year")
    print("m3")
    sanity(m3)
    dd$m3_res <- residuals(m3)

    # 4. logistic oxygen
    # m4 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + temp_sc + temp_sq + logistic(oxy_sc),
    #              data = dd,
    #              mesh = mesh,
    #              family = delta_gamma(link1 = "logit", link2 = "log"),
    #              spatiotemporal = "off",
    #              spatial = "on",
    #              time = "year",
    #              control = sdmTMBcontrol(newton_loops = 2))
    # print("m4")
    # sanity(m4)
    # dd$m4_res <- residuals(m4)

    # 5. linear metabolic index
    m5 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + phi_sc,
                 data = dd,
                 mesh = mesh,
                 family = tweedie(link = "log"),
                 spatiotemporal = "IID",
                 spatial = "off",
                 spatial_varying = ~0 + quarter_f,
                 time = "year")
    print("m5")
    sanity(m5)
    dd$m5_res <- residuals(m5)

    # 6. breakpoint metabolic index!
    m6 <- sdmTMB(density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + breakpt(phi_sc),
                 data = dd,
                 mesh = mesh,
                 family = tweedie(link = "log"),
                 spatiotemporal = "IID",
                 spatial = "off",
                 spatial_varying = ~0 + quarter_f,
                 time = "year")
    print("m6")
    sanity(m6)
    dd$m6_res <- residuals(m6)
    
    # logistic metabolic index!
    # m7 <- update(m1, density ~ 0 + year_f + quarter_f + sal_sc + depth_sc + depth_sq + logistic(phi_sc))
    # print("m7")
    # sanity(m7)
    # dd$m7_res <- residuals(m7)
    
    # Plot residuals
    p1 <- dd %>%
      dplyr::select(m0_res, m1_res, m2_res, m3_res, m5_res, m6_res) %>%
      pivot_longer(everything()) %>%
      ggplot(aes(sample = value)) +
      stat_qq(size = 0.75, shape = 21, fill = NA) +
      facet_wrap(~name) +
      stat_qq_line() +
      labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
      theme(aspect.ratio = 1)
    
    print(p1)

    # Save AIC
    data_list_aic[[i]] <- AIC(m0, m1, m2, m3, m5, m6) %>%
      tibble::rownames_to_column("model") %>%
      mutate(group = i)
    
}
filter: removed 54,050 rows (81%), 12,336 rows remaining
Loading required namespace: INLA
[1] "cod_adult"
[1] "m0"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m1"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m2"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m3"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m5"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m6"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
pivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 12336x6, now 74016x2]

mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 54,050 rows (81%), 12,336 rows remaining

[1] "cod_juvenile"
[1] "m0"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m1"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m2"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m3"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m5"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m6"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
pivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 12336x6, now 74016x2]

mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 54,802 rows (83%), 11,584 rows remaining

[1] "plaice_adult"
[1] "m0"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m1"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m2"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m3"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m5"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m6"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
pivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 11584x6, now 69504x2]

mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 54,802 rows (83%), 11,584 rows remaining

[1] "plaice_juvenile"
[1] "m0"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m1"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m2"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m3"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m5"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m6"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
pivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 11584x6, now 69504x2]

mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 57,114 rows (86%), 9,272 rows remaining

[1] "flounder_adult"
[1] "m0"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m1"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m2"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m3"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m5"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m6"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
pivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 9272x6, now 55632x2]

mutate: new variable 'group' (character) with one unique value and 0% NA
filter: removed 57,112 rows (86%), 9,274 rows remaining

[1] "flounder_juvenile"
[1] "m0"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m1"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m2"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m3"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m5"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
[1] "m6"
✔ Non-linear minimizer suggests successful convergence
✔ Hessian matrix is positive definite
✔ No extreme or very small eigenvalues detected
✔ No gradients with respect to fixed effects are >= 0.001
✔ No fixed-effect standard errors are NA
✔ No standard errors look unreasonably large
✔ No sigma parameters are < 0.01
✔ No sigma parameters are > 100
✔ Range parameter doesn't look unreasonably large
pivot_longer: reorganized (m0_res, m1_res, m2_res, m3_res, m5_res, …) into (name, value) [was 9274x6, now 55644x2]

mutate: new variable 'group' (character) with one unique value and 0% NA

# Save aic as data frames
data_aic <- bind_rows(data_list_aic)

write_csv(data_aic, paste0(home, "/output/data_aic_01.csv"))

Plot meshes in a single plot using patchwork

g1 <- d %>% filter(group == unique(d$group)[1])
filter: removed 54,050 rows (81%), 12,336 rows remaining
mesh_g1 <- make_mesh(g1, xy_cols = c("X", "Y"), cutoff = 15)

g2 <- d %>% filter(group == unique(d$group)[5])
filter: removed 57,114 rows (86%), 9,272 rows remaining
mesh_g2 <- make_mesh(g2, xy_cols = c("X", "Y"), cutoff = 15)

g3 <- d %>% filter(group == unique(d$group)[3])
filter: removed 54,802 rows (83%), 11,584 rows remaining
mesh_g3 <- make_mesh(g3, xy_cols = c("X", "Y"), cutoff = 15)

g4 <- d %>% filter(group == unique(d$group)[2])
filter: removed 54,050 rows (81%), 12,336 rows remaining
mesh_g4 <- make_mesh(g4, xy_cols = c("X", "Y"), cutoff = 15)

g5 <- d %>% filter(group == unique(d$group)[6])
filter: removed 57,112 rows (86%), 9,274 rows remaining
mesh_g5 <- make_mesh(g5, xy_cols = c("X", "Y"), cutoff = 15)

g6 <- d %>% filter(group == unique(d$group)[4])
filter: removed 54,802 rows (83%), 11,584 rows remaining
mesh_g6 <- make_mesh(g6, xy_cols = c("X", "Y"), cutoff = 15)


p1 <- ggplot() +
  inlabru::gg(mesh_g1$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = g1, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g1$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) + 
  labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g1$group[1], "_", " "))) +
  theme_sleek(base_size = 8) +
  theme(axis.title.x = element_blank(),
        plot.margin = unit(c(0,0,0,0), "cm"))
  
p2 <- ggplot() +
  inlabru::gg(mesh_g2$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = g2, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g2$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) + 
  labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g2$group[], "_", " "))) +
  theme_sleek(base_size = 8) +
  theme(axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        plot.margin = unit(c(0,0,0,0), "cm"))

p3 <- ggplot() +
  inlabru::gg(mesh_g3$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = g3, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g3$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) + 
  labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g3$group[], "_", " "))) +
  theme_sleek(base_size = 8) +
  theme(axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.title.y = element_blank(),
        plot.margin = unit(c(0,0,0,0), "cm"))

p4 <- ggplot() +
  inlabru::gg(mesh_g4$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = g4, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g4$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) + 
  labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g4$group[], "_", " "))) +
  theme_sleek(base_size = 8) +
  theme(axis.title.x = element_blank(),
        plot.margin = unit(c(0,0,0,0), "cm"))

p5 <- ggplot() +
  inlabru::gg(mesh_g5$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = g5, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g5$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) + 
  labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g5$group[], "_", " "))) +
  theme_sleek(base_size = 8) +
  theme(axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = unit(c(0,0,0,0), "cm"))

p6 <- ggplot() +
  inlabru::gg(mesh_g6$mesh) +
  coord_fixed() +
  geom_point(aes(X, Y), data = g6, alpha = 0.2, size = 0.5) +
  annotate("text", -Inf, Inf, label = paste("n knots = ", mesh_g6$mesh$n), hjust = -0.1, vjust = 2, size = 2.5) + 
  labs(x = "Easting (km)", y = "Northing (km)", title = str_to_sentence(str_replace(g6$group[], "_", " "))) +
  theme_sleek(base_size = 8) +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = unit(c(0,0,0,0), "cm"))

p1 + p2 + p3 + p4 + p5 + p6

ggsave(paste0(home, "/figures/supp/mesh.pdf"), width = 17, height = 11, units = "cm", device = cairo_pdf)